Phase field model of solid-liquid and liquid-liquid phase transitions 
in flow and elastic fields in one-component systems 
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Abstract. - We construct a phase field model including hydrodynamics and elasticity in one- 
component systems. It can be used to investigate solid-liquid and liquid-liquid phase transitions. 
Upon first-order phase transition, a velocity field is induced around interfaces in the presence of 
a density difference between the two phases even without applied shear flow. As applications, we 
present simulation results on two cases of melting, where a solid domain is placed on a heated wall 
in one case and is suspended in a warmer liquid under shear flow in the other. We find that the 
solid domain moves or rotates as a whole due to elasticity, releasing latent heat. We also examine 
the liquid-liquid phase transition of a highly viscous domain into a less viscous liquid on a heated 
wall, where an inhomogeneous velocity field is induced within a projected part of the domain. In 
these phase transitions, the interface temperature is nearly equal to the coexisting temperature 
Tcx{p) away from the heated wall in the presence of heat flow in the surrounding liquid. 



Introduction. — In solid-liquid phase transitions, 
various nonequilibrium patterns have been observed [Tl[2] . 
To reproduce such patterns numerically, phase field mod- 
els have been used extensively [3l|4], where a space-time 
. dependent phase field 0(r, t) takes different values in solid 
I and liquid varying smoothly across diffuse interfaces. A 
. merit of this approach is that any surface boundary condi- 
■ tions need not be imposed explicitly in simulations. Most 
theories of crystal growth have assumed that the dynam- 
ics of first-order phase transition is governed by diffu- 
. sion of heat and/or composition. Though some attempts 
' have been made to include the velocity field or convection 
into theory [SHE], understanding of hydrodynamic effects 
during solidification or melting still remains inadequate. 
Moreover, phase field calculations were performed on the 
surface instability in epitaxial film growth [9l|T0], where 
elastic effects are crucial. It is worth noting that various 
phase field models have been used to investigate phase 
ordering in solid-solid phase transitions pT] . 

We also mention phase transitions in one component flu- 
ids between two Hquid phases with different microscopic 
structures JSHS]- In particular, Tanaka's group ^ per- 
formed experiments of the phase ordering dynamics at 
Hquid-liquid phase transitions. To interpret their data, 
they introduced a nonconserved order parameter repre- 



senting microscopic structural order. From our viewpoint, 
we should develop a phase field model of liquid-liquid 
phase transitions, where the structural order parameter, 
denoted by the same notation (/), is coupled to the hydro- 
dynamic variables. 

One of the present authors has developed the dynamic 
van der Waals theory for one-component fluids p^[T5] . It 
is a phase field model of fluids based on the van der Waals 
theory. It can treat evaporation and condensation with 
inhomogeneous temperature accounting for latent heat. 
The aim of this letter is to present a phase field model 
including both hydrodynamics and elasticity on the ba- 
sis of well-defined thermodynamics. Our model should 
thus be applicable to solid-liquid and liquid-liquid phase 
transitions. We will treat one-component systems only for 
simplicity. No gravity will be assumed. 

Phase-field model. — We discuss thermodynamics 
and dynamics for the phase field (j) and the hydrodynamic 
variables. In the bulk region (j) is equal to in liquid (or 
liquid phase I) and to 1 in solid (or liquid phase II). 

At the starting point, we introduce an entropy density 
including a gradient contribution as 



5(n,e,(/.)-ic|V0|2, 



(1) 
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where 6'(n, e, 0) is a function of the number density n, 
the internal energy density e, and (j). In this work C is 
a positive constant (which may depend on (j) more gener- 
ahy). The temperature T and the chemical potential ji 
are defined by 1/T = dS/de and /i/T = -dS/dn. The 
derivative with respect to (j) is written as dS/d(j) = —T/T. 
The differential form of S reads 



dS = {de - iidn - Tdcj)) /T. 



(2) 



Neglecting the gradient energy density jTl] , we assume the 
total energy density in the form, 



ex = e + p\v\^/2 + G((/))$(e2, eg), 



(3) 



where p = mn is the mass density with m being the molec- 
ular mass, V is the velocity field, and the last term is the 
elastic energy density. The shear modulus G{(j)) is zero for 
= (in liquid) and is a positive constant for = 1 (in 
solid) for solid-liquid transitions, while G = for liquid- 
liquid transitions. Here we suppose two dimensions, where 
anisotropic elastic strains are 62 and 63. For small elas- 
tic deformations in solid, we have 62 = VxUx — ^yUy^ 
63 = ^xUy + VyUx^ and ^ = (62 + 63)72 in terms of the 
displacement field u = {ux^Uy). Hereafter V^; = d/dx and 
Vy = d/dy. In the simulation in this work, the amplitudes 
of 62 and 63 remain very small (< 10 ~^). However, the 
linear elasticity does not hold for large strains and a sim- 
ple form of ^ applicable in the nonlinear regime is given 
by the periodic function ^16^ 



1 

6^ 



3 — cos(27re2) — cos(27re+) — cos(27re_) 



(4) 



where e± = (>/3e3 ± e2)/2. If we set 63 = ecosx and 62 = 
esinx, we have ^ = jl-i^^e^ j^-i^^e^ cos(6x)/720+- • •. 
If we rotate the reference frame by angle ^, X is shifted by 
2d. Thus <I> is highly isotropic for small e (less than 0.5). 

We follow the principle of nonnegative entropy produc- 
tion to set up the dynamic equations. We introduce the 
generalized thermodynamic force associated with ^ as 



r = r - Tcv^ 



(5) 



where G' = dG/d4>- The dynamic equation of </> is then 



dt 



Ar, 



(6) 



where A is the kinetic coefficient. The hydrodynamic 
equations are of the usual forms. 



^ = -V • {pvv+ n 
ot 



^ = -V • [ctv + {11- a) -v- AVT]. 



(7) 
(8) 
(9) 



The reversible stress tensor 11= 
ent and elastic parts as 



{I^ij} contains the gradi- 



H,, = {p-G^--TC\V^\^)8iJ^TCV^^VJ^-GeiJ. (10) 



where p = nji — e 
sions, the tensor e^^ 

and Cnr-n, ^IIT, 



ST is the pressure. In two dimen- 



is given by by 



^yy 



d^/de2 



d^/des' The viscous stress tensor 



^= Wij} = {vC^i'^j + Vj^i) + {rjB - 2r]/d)5ijV • v} is 
written in terms of the shear viscosity 77 and the bulk vis- 
cosity rjB (in d dimensions). The strains 62 and 63 obey 



de2 
~dt 
des 

dt 



-V • Ve2 + VxVx - "^yVy, 
-V • Ve3 + \/xVy + "^yVx- 



(11) 

(12) 



With these equations, the entropy density S obeys 



-V- 



Sv - AtCV4) - ^VT 

- e« + e<t,)/T, 



(13) 



where e„, eg, and are the heat production rates, 



A 



-(VT)2, e„ = ^ai,Vit;,-, = Af^. (14) 



Model entropy. — We suppose a reference equilib- 
rium state at T = To and p = Po = Pcx{To)^ where liquid 
and solid coexist macroscopically and the chemical poten- 
tial takes a common value p = fiQ = /icx(^o)- The quanti- 
ties in the reference liquid state will be denoted with the 
subscript 0^, while those in the reference solid state with 
the subscript Os. The number and energy densities in the 
reference liquid (solid) are written as no£ and eo£ (nos and 
eos) , respectively. The number and energy density devia- 
tions are defined as those from the reference liquid values 
as (5n = n — no£ and 6e = e — eo£. 

We propose to use a simple expression for the entropy 
density 5* in Eq.(2.2). It contains terms up to second or- 
ders in 5n and Se as 



S = So£ 



Se — poSn 
To 



2ToK° 



W{<P). (15) 



Here So£^ Cy, and K!^ are the entropy density, the 
const ant- volume heat capacity per unit volume, and the 
isothermal compressibility in the reference liquid state, re- 
spectively. We define dimensionless variables r and ( by 

r = {5e-f3oSn)/ToC^-aiO{^), (16) 
C = 5n/nQ£- a20{(j)), (17) 

where /3o = {de/dn)T is the derivative in the reference 
liquid state. We define VK(^) in S and 0{(j)) in r and C as 



w{4>) = \a^\i 



■4>f, e{4>) = 4>\i-'^4>), 



(18) 
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where A is a positive constant. Then ^(0) =0 and ^(1) = 
1. The coefficients ai and a2 in r and are related to 
the latent heat and the density difference between liquid 
and solid, as will be shown below. For liquid with ^ = 0, 
Eq.(15) is a well-known expansion form up to second order 
deviations ^Tj. Adding the terms proportional to O^cj)) in 
r and we use it even for solid. 

As functions of r, and (/), T, /i, and T are written as 

T = To/(1-t), (19) 
M = [^Jio-PoT + C|nuK^]/{l-T), (20) 

r = Tj^[W{<P)-nMhem, (21) 

where h is the linear combination of r and ^ defined by 

noih = (aiC^)r + {a2/kBToK^)C. (22) 
The pressure p = nfi-\- TS — e is then of the form, 



P-Po 



c 



-W^noihO, (23) 



where = —{dn/dT)p/n is the thermal expansion coef- 
ficient in the reference liquid state. 

We consider equilibrium two-phase coexistence without 
anisotropic strains. From Eqs.(19) and (20) r and ( take 
common values in the two phases, while Eqs.(21) and (23) 
yield h = 0. The liquid phase ^ = is preferred for < 
and the solid phase = 1 for > 0. Thus h represents 
the distance from the coexistence line in the phase dia- 
gram. From Eq.(17) the density difference between the 
two phases is given by 



An : 



Hi = a2noi. 



(24) 



The latent heat per particle is q = T{s£ — 5^), where s = 
S/n. Its value in the reference state (r = 0) is written as 



qo = BiToa2/noi{l + a2). 



(25) 



where Bi = a^/K^ — Cyai/a2 is equal to the deriva- 
tive {dp/dT)cx along the coexistence line for r = 0. 
If h — 0, we conffim the Clapeyron-Clausius relation 
q = T{dp/dT)cx^n/n£ns for any r. On the other hand, 
the equation for a planar interface (j) = (j)int{z) becomes 
VK(^) = C\d(j)/dz\^^ which is solved to give (j)int{z) = 
1/(1 + e^/^). The interface thickness is ^ = (C/A)^/^ and 
the surface tension is 7 = 

Numerical Method. — In two examples to follow, we 
integrated Eqs. (7), (8), and (11)- (13) in two dimensions 
on a 800 x 400 lattice. We use the entropy equation (13) 
instead of the energy equation (9) to suppress the so-called 
parasitic flow [15], which is an artificial flow around an 
interface pTtfTS] . The simulation mesh length is Ax = ^ = 
(C/A)^/^ Here we set A = OAnoikB, so C = OAnoikB^^- 
The horizontal and vertical lengths of the system are then 
L = 800^ and H = 400^, respectively. We imposed the 




^0.008^/to solid: 32 = 0.2 



t = 2000 to 




0.008^ /to solid: 32 = - 0.2 



Fig. 1: Solid domain at four times on a heated wall for a2 = 0.2 
(upper four plates) and for a2 = —0.2 (lower four plates) in 
the region 1/2 < x/L < 7/8 and < y/H < 3/4. Arrows 
represent the velocity field. Those below the panels correspond 
to 0.008^/to. Melting mostly takes place close to the heated 
wall, where the velocity is from solid to liquid for a2 = 0.2 and 
from liquid to solid for a2 = —0.2. 



periodic boundary condition along the horizontal x axis 
and the no-slip condition t; = at = and H. In 
addition, neglecting the surface entropy and energy, we 
set d(j)/dz = at y = and H. 

The thermodynamic quantities are given by = 
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— 0.008^ /to fluid: Mv = 50 32 = 0.2 

Fig. 2: Droplet at liquid- liquid transition with C = at t/to = 
2000 and 5000 on a heated wall in the region 1/2 < x/L < 7/8 
and < yjH < 3/4. The viscosity ratio My is 5 (top) and 
50 (bottom). A typical velocity is 0.02^/to for My = 5 and 
0.008^/to for My = 50. A circulating velocity field is marked. 
See Fig. 5 for the velocity along the interface. 



t = 2000to t/Tq 



t = 5000tn T/To 




0.8 0.7 6 5 0.6 
x/L x/L 

solid: 32 = 0.2 



Fig. 3: Temperature in the region 0.5 < x/L < 0.8 and < 
y/H < 1 for the solid case with a2 = 0.2. See the corresponding 
domain shapes in Fig.l. At t = 2000to (left), heat from the wall 
is used to melt the solid near the contact point. At t = 5000to 
(right), the domain is detached from the wall and its interface 
temperature becomes a constant. 




0.4 
x/L 



y/H 

solid: 32 = 0.2 



Fig. 4: h in Eq.(22) of a melting solid along the interface, which 
is multiplied by 100 at t = 2000to (left) and by 1000 at t = 
5000to (right). It represents the distance from the coexistence 
curve, where /i > in the body part (blue) and /i < in the 
constricted part (red). 



I.IC^ = 9A7noikB, = OMA/noikBTo, and = 
0.20/To, which are consistent with the thermodynamic 
identity Ta^ = {Cp - Cy)KT of fluids [TlJ. The shear 
modulus is G = 10niokBTo<p'^ for solid-liquid transitions. 
By setting ai = —0.27, we examine two cases of a2 = 
An/no£ = 0.2 and —0.2. The latent heat qo in Eq.(25) is 
2.5kBTo for as = 0.2 and IMbTq for as = -0.2. 

The viscosities are given by r]{(j)) = r]B{<i)) = + {Vs — 
r]i)(p'^. Let = rji/mriQi be the liquid kinematic viscosity. 
We measure space and time in units of ^ and 



to = 0.4^V^o- 



(26) 



The viscosity ratio My = r]s/Vi is taken to be 10 or 50. 
The simulation time mesh At is O.OOSto for My = 10 and 
0.0025^0 for My = 50. The thermal conductivity is A = 
1.63i'oCp{l + 3^^), so the Prandtl number is 1/1.63 in 
liquid. The kinetic coefficient for ^ is A = 0.16/no^/cB^o^o- 
With these expressions, the reference temperature Tq and 
density no£ do not appear in the scaled dynamic equations. 

A semispheric domain on a heated substrate. — 

First, we placed a solid (type II liquid) semisphere with 
radius R = 200^ on the substrate in liquid (type I liquid), 
where cj) = 1 inside the domain, ^ = outside it, and 
T = C = 62 = 63 = in the whole cell. Thus T = Tq and 
Sn = noia20{(j)) were imposed. The boundary tempera- 
tures at z = and H were fixed at Tq. Then small re- 
laxations followed near the interface in short times (~ to). 
After an equilibration time of lOOto, we raised the bottom 
temperature to 1.05To. We set t = at this heating. 

In Fig.l, we show the domain shapes and the surround- 
ing velocity at four times. The solid density is higher or 
lower than the liquid density depending on the sign of a2. 
The phase change takes place mostly near the heated wall. 
For a2 = 0.2 the flow is from the solid to the liquid near the 
bottom in the upper plates, while for a2 = —0.2 the flow is 
in the reverse direction in the lower plates. The solid veloc- 
ity is very small and is nearly uniform within the domain. 
The strains 62 and 63 remain of order 10~^ and are well 
in the linear elasticity regime. Nevertheless, the resultant 
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Fig. 5: Liquid velocity Vi, interface velocity Vint, and melting 
flux J = ni{vi—Vint) in the normal direction along the interface 
in units of ^/to or no^^/to, where t = 2000to and a2 — 0.2. 
They are for liquid-liquid transitions with Mv = 5 (top) and 
50 (middle) and for a solid- liquid transition (bottom). Melting 
occurs mostly in the constricted part close to the heated wall. 



shear stress ('^ lO^^no^/cBTb) can realize the solid body 
motion. In addition, the pressure p in the liquid region 
gradually increases for a2 > and decreases for a2 < due 
to the density difference. In our examples the deviation 
p — pQiii the liquid is of order iClSno^/cBTb for a2 = ±0.2 
at t = IOHq. (Here the resultant adiabatic temperature 
change is very small since {dT/dp)s ~ 0.03/no^/cB-) 

Figure 2 gives the profiles of a droplet in a less vis- 
cous liquid at liquid-liquid transitions with G = 0. The 
velocity field decreases with increasing the viscosity ratio 
My = r]s/r]i, but it still remains noticeable in the projected 
part of the domain even for My = 50. Here an interface 
motion induces a fluid motion with a small velocity gradi- 
ent within a highly viscous droplet. However, the domain 
shapes in the two cases in Figs.l and 2 are very similar. 



In Fig. 3, the temperature is displayed for the solid case 
with a2 = 0.2 before and after detachment of the domain. 
Similar profiles were also obtained for the case of liquid- 
liquid transitions. A steep gradient can be seen near the 
heated wall, but the interface temperature is nearly flat at 
the melting temperature T^^ip) far from the wall. Thus 
the interface is divided into the body part far from the 
wall and the constricted part close to the wall. In this 
case the thermal diffusion length {DrtY^'^ ~ (V^o)^^^^ is 
still shorter than the cell height H = 400^, resulting in a 
small temperature gradient in the upper region. In Fig. 4, 
we show that h in Eq.(22) is very small on the interface. 
Here h ^ -2.b[r - {dT/dp)c^{p - po)/Tc] away from the 
interface (or for = or 1 in Eq.(23)). In the constricted 
part, the gradient of h is —10% of that of r. 

In Fig. 5, we plot the liquid velocity V£ = u • the 
interface velocity Vint = • ^int-, and the melting flux J = 
niivi — Vint) through the interface, where is the liquid 
density. These are the quantities close to the interface in 
the normal direction ly. Indeed, melting mostly occurs in 
the constricted part for large My and for the sohd case. 
We also notice that \vint\ is considerably larger than l^^l, 
which is obviously due to the small size of a2. 

A solid domain in a sheared warmer liquid. — 

Next, at t = 0, we placed a 240x240 solid square with 
a2 = 0.2 in a liquid. The temperature was initially 0.95To 
in the solid and I.ITq in the liquid. The top and bottom 
plates were insulating or dT/dz = 0. By moving them we 
applied a shear flow with rate 7 = 2 x 10~Hq^ for t > 0. 
The liquid should be cooled upon melting. 

Figure 6 displays the solid shapes at t = lO^to = 0-2/7 
and lOHo = 2/7 in the middle region 1/4 < x/L < 3/4 
and < y/H < 1, where the solid is rotating and melt- 
ing. The final solid area in Fig. 6 is 70% of the initial 
area. Figure 7 gives the temperature profile at t = IOHq. 
It demonstrates a considerable cooling of the liquid par- 
ticularly in the regions where the melted liquid has been 
convected. The temperature is homogeneous with h = 
within the solid domain. The average temperature (T) = 
J^dx dyT{r,t)/LH is 1.078To at t = 0, 1.076To at 
t = lOHo, and 1.063To at t = 10%. To explain the cool- 
ing by O.OI5T0 at t = lOHo, we multiply the number of 
melted particles (~ 4000) by the latent heat per particle 
qo in Eq.(25) and divide it by the total specific heat LHC^ 
to obtain a temperature decrease ~ 0.015To in accord with 
the data of (T). 

Summary. — In our theory, the entropy depends on 
the phase field (/) and contains a gradient part (oc |V^p), 
while the total energy consists of the usual internal energy, 
the kinetic energy, and the elastic energy. The elasticity is 
introduced using strain fields 62 and 63 in two dimensions, 
which represent anisotropic elastic deformations. Starting 
with these quantities and using the principle of nonneg- 
ative entropy production, we have constructed dynamic 
equations for the phase field, the hydrodynamic variables. 
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Fig. 6: 



A square-shaped solid domain in , 

-4,-1 



, warmer liquid under 



shear flow with 7 = 2 x 10~^t^^ at = 0.2 (top) and 2 (bot- 
tom). It gradually melts and rotates as a solid body. Melted 
particles are convected to cool the surrounding liquid. 



lOOOOto 




Fig. 7: Temperature profile around a melting solid in shear flow 
at 7t = 2, corresponding to the lower panel of Fig. 6. It is lowest 
in the directions making angles of 7r/4 and 57r/4 with respect 
to the X axis due to convection of cooler melted particles. 



and the strains. They can describe phase transition dy- 
namics accounting for the hydrodynamic and elastic ef- 
fects. If the strains are absent or the shear modulus 
vanishes, we obtain the dynamic equations applicable to 
liquid-liquid phase transitions. We have also proposed 
to use the entropy density containing second-order devi- 
ations of the number and internal energy densities from 
a reference two-phase state. We have solved the dynamic 
equations to examine two melting phenomena in two di- 
mensions. In these cases small strains of order 10 ~^ have 
produced rigid body motions of a solid domain. 

There can be a number of problems to be studied in our 
scheme such as dendrite formation in flow, spinodal de- 
composition, epitaxial growth, and recrystallization. In 
particular, it is of great interest to investigate kinetics 
with large strains or even with dislocations. In liquid- 
liquid transitions, we should further investigate the effects 
of latent heat and shear flow in phase ordering. 
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